Torsional vibration analysis of shaft with multi inertias

An analytical method is proposed to investigate the torsional vibration of the uniform circular shaft with multiple concentrated inertias. The governing equation is established based on the Hamiltonian principle and verified by the dynamical method. The theoretical solutions of frequencies and mode shapes under different boundary conditions are obtained using the separation variable method and integral transformation. The effectiveness of the proposed method is verified by comparison with existing literature. Considering the change of the magnitudes/positions/number of concentrated inertias, and different boundary conditions, the natural frequencies and mode shapes are discussed. Several general rules are obtained. Moreover, some interesting phenomena have been found and explained. The analytical method has applications in the design of shafting with multiple concentrated inertias and the reliability checking of the “approximate” solutions.

gave the exact expressions for the torsional frequencies and mode shapes of the generally constrained shafts and piping with two disks. Reference 19 studied the torsional vibration problem of a beam containing lumped elements under cantilever boundary conditions and fixed boundary conditions. As for the "exact" solutions of the torsional vibration of uniform shafts carrying multi concentrated elements, the information concerned is rare. Only Chen 20 presented an "exact" solution for free torsional vibration of a uniform circular shaft carrying multiple concentrated elements adopting the numerical assembly method (NAM). But NAM [21][22][23] does not give the "exact" expression of mode shapes, and the method needs to establish complex coefficient matrices, which is still quite complicated. In order to remove the barriers, an analytical method (AM) based on the Hamilton principle 6,24 , variable separation method, and integral transformation is used to perform the torsional vibration analysis of a uniform circular shaft carrying multiple concentrated elements (rotary inertias) with arbitrary magnitudes and locations. The paper is organized as follows: In "Governing equations" section, the simplified model, the establishment and verification of the governing equation are introduced. "Solution of the torsional vibration equation" section gives the solution of the governing equation with different boundary conditions, including the analytical expressions of mode shapes natural frequencies. "Validation of the proposed method" section shows the validation of the proposed method compared with the existing literature. The applications of AM are presented in "Application of AM" section, and several special phenomena are discussed in "Discussion of some phenomena" section. Finally, in "Conclusion" section, the concluding remarks and applications are discussed.

Governing equations
Model of the problem. The physical model is shown in Fig. 1 by assuming the cross-section remains flat during torsion and the lateral deformation of the shaft is neglected. Therefore, only the torsional deformation of the shaft is considered. Governing equations established by the Hamiltonian principle. The angular displacement of each cross-section is θ = θ(x, t) and γ = r∂θ ∂x is the shear strain at each point in the cross-section. Hence, considering the constitutive relation, the shear stress at each point in the cross-section is τ = Gγ.
The deformation energy of the shaft can be expressed as where I P (x) represents the moment of inertia of the cross-section. Considering the problem of uniform cross-section shaft I P (x) is constant which degenerates into I P .
The kinetic energy of the shaft obeys The work of torque has the form of  (1)-(4) into Eq. (5) and considering δθ| t=t 0 = 0 and δθ| t=t 1 = 0 , energy variation is shown as Equation (6) is suitable for arbitrary t 0 , t 1 . Then, δθ, δθ| x=l , δθ| x=0 are arbitrary. Thus, the vibration differential equation can be expressed as The boundary conditions are written as When the concentrated inertia is not taken into account, the equation is simplified into the differential equation of free vibration of the continuous beam Considering free-boundary conditions (free-free, FF), Eq. (11) can be expressed as When the torque is zero and the boundary are fixed (pinned-pinned, PP), Supposing M 0 = M l = 0 , the mixed boundary condition (free-pinned, FP) is given as The verification of governing equations. To verify the correctness of the vibration equation established in "Governing equations established by the Hamiltonian principle" section, it is analyzed from the perspective of torque balance. Figure 2 shows a beam element without concentrated elements (rotary inertias).
∂θ ∂x x=0,l = 0.  www.nature.com/scientificreports/ The length of the beam element is dx , the torque on its left is M and on the right side is M + ∂M ∂xdx . Figure 3 shows a beam element with the i-th concentrated rotary inertia located at x = x m,i , and has the same geometric parameters as Fig. 2.
When the micro-segment shown in Fig. 2 does not contain the concentrated rotary inertia, the equilibrium equation can be obtained When the micro-section contains the i-th concentrated inertia, Eq. (12) becomes into Supposing multiple concentrated inertias, the equation expression obeys where Substituting Eq. (17) into Eq. (16), the differential equation of vibration is the same as Eq. (7). In the same way, the boundary conditions can be obtained, which are coincident with Eqs. (8), (10)- (13). Therefore, the theory applied in this paper is credible, and the equations established are also available.

Solution of the torsional vibration equation
.. .

Validation of the proposed method
In order to validate the proposed method, three cases were studied and compared with the results published in Refs. 6,16,20 . Case one. The cantilevered shaft carrying single rotary inertia, located at x m,1 /l = 1.0 or 0.5 respectively, were studied first (Fig. 4). The dimensions and physical properties of the circular shaft studied in case one are the same as those in Ref. 20 (29). Then non-dimensional frequency coefficients can be obtained. When ξ 1 = 1 , the Eigen equation is given as Eq. (31), which is the same as it in Ref. 6 .
When ξ 1 = 0.5 , the Eigen equation is given as  16 and Chen 20 , the lowest five non-dimensional frequency coefficients, β j (j = 1-5) were shown in Table 1. Table 1 shows that the results of AM, Gorman 16 , and Chen 20 are in good agreement, which means that the proposed method is accurate and effective for a simple model. Case two. Further, to verify the applicability of the proposed method in complex models, case two was adopted. Hereby, a shaft carrying five inertias, as shown in Fig Table 2, and the corresponding mode shapes of twisting angles were shown in Figs. 6 and 7.
The differences between AM and � Reference (� FEM , � Chen ) shown in the parentheses ( ) of Table 2 were calculated with the formula: ε j = (� Reference, j − � AM, j )/� AM,j , where FEM,j and Chen,j denote the j-th natural frequencies of the shafts carrying ''five'' torsional springs obtained from AM, Chen 20 and FEM, respectively. From Table 2 one finds that values of AM,j agreement with FEM,j and Chen,j (for PP and FP shafts), hence the accuracy of the AM is good. It also can be found that values of ε j between Chen,j and AM are all smaller than those between FEM,j and AM . This is because the FEM model is approximate, and the calculation accuracy of natural frequency can be improved by improving the degree of freedom of the FEM model. Figures 6 and 7 show the lowest five mode shapes of twisting angles for the shaft carrying five concentrated elements in the support conditions: PP and PF, obtained from AM and Chen 20 . Mode shapes' variation characteristics with x/l , obtained by AM, are almost identical to the ones calculated by Chen 20 . However, for the same order mode, the mode shape obtained via AM is slightly different from Chen's at the points where rotary inertias  www.nature.com/scientificreports/ exist, which should be blamed on the application of the conclusion from Chen 20 and the way of mapping. In Ref. 20 it was pointed out that mode shapes of the shafts carrying five inertias are almost coincident with the ones for the shaft without carrying any concentrated elements and the former was not shown. To compare the results of the mode shapes, the aforementioned conclusion was directly used in this paper. The mapping way of mode shape was combing limited discrete points with trend lines. As a result, if the number of discrete points is small, a similar conclusion to Ref. 20 can be obtained via AM, as shown in Figs. 8 and 9. Once the number of discrete points is increased, the aforementioned conclusion is no longer accurate (shown in Figs. 6, 7). In fact, concentrated elements have a significant influence on mode shapes, which can be found in many pieces of literatures 7,21,[25][26][27] .
To verify the accuracy of the mode shapes obtained via the presented method, FEM was employed (solved by Ansys 14.0 https:// www. ansys. com/ zh-cn/ produ cts/ struc tures/ ansys-mecha nical). The comparison of the lowest five mode shapes obtained by FEM with 51-element uniform discrete and AM was given in Figs. 10 and 11. While modeling by FEM, a simplified model of a uniform discrete section should include a sufficient number of elements with higher-order interpolation functions. Here, beam188 with quadric shape functions for 3-D (3-node) line element was applied. From Figs. 10 and 11, it can be seen that AM results are in good agreement with FEM results, indicating that the function of the mode shape via the proposed method is accurate.
The findings in Table 2 and Figs. 6,7,8,9,10 and 11 indicate that the proposed method is accurate and effective for complex models. The results also show that concentrated elements have a significant influence on the modal shape. This phenomenon has also been found in Ref. 7,21,[25][26][27] , which is another evidence of the accuracy of AM. The reason for the phenomenon is given in "Explanation of Phenomenon 5" section.
Case three. To determine the engineering applicability of the method, case three was employed. A free boundary shaft with one moving concentrated inertia, as shown in Fig. 12 Fig. 13.
From Fig. 13, it is seen that the lowest six natural frequencies of the FF shaft carrying one concentrated inertia analyzed by AM are the same as those obtained by LMM, which suggests that the presented method is applicable in practice.
According to the aforementioned comparisons, it is believed that the AM is an effective way for the title problem.          The concentrated inertia located at a certain point. To further investigate the phenomenon in "Position influence of the concentrated inertia" and "Value influence of the concentrated inertia" sections, the model in case three continued to be studied. In this case, J m has a certain position ( x m = 0.5l, 0.25l, 0.279l, and 0.679l ) and a changing radius R(0 ≤ R ≤ 0.2m) . Natural frequencies (1st-6th) were given in Fig. 15, and mode shapes were given in Figs. 16 and 17.    Fig. 15, when J m is located at x m = l/n 1 (n 1 = 2, 4, 6, . . .) , the n 1 (2q + 1)/2 (q = 0, 1, 2, . . .)-th natural frequency does not change with R(J m ) , the other reduce with R(J m ) , while each frequency trends to a certain value and the n 1 (2q + 1)/2 + 1-th natural frequency equals to the n 1 (2q + 1)/2-th natural frequency if R(J m ) → ∞ ; when J m is not located at x m = l/n 1 (n 1 = 2, 4, 6, . . .) , each natural frequency reduces with R(J m ) and each frequency trends to a certain value if R(J m ) → ∞ . The mentioned phenomena also influence the mode shapes of the shaft. Therefore, Figs. 16 and 17 were given for further discussion.
In Fig. 16, mode shapes of the lowest six natural frequencies were given when J m is located at x m = l/2 . It can be seen that the n 1 (2q + 1)/2 (n 1 = 2, q = 0, 1, 2, . . .)-th order mode shape did not change with R , the other mode shapes change evidently and the range of those change mode shapes are beyond the interval [− 1, 1] 28-30 . From Fig. 17 the same conclusions can be obtained.

Number influence of the concentrated inertias. Number influence of concentrated inertias on
the natural frequency was investigated in this case ( Table 3). The results were given in Fig. 18. N m = 1−4 , 0 ≤ R ≤ 0.4 m , f j j = 1, 2, 3, . . . : the j-th natural frequency and f j O : the j-th natural frequency, when R = 0. Figure 18 shows that each natural frequency reduces when values of the concentrated inertias increase or when the number of concentrated inertias increases. These findings indicate that the natural frequency can be controlled by changing the magnitude and number of inertia.
The influence of boundary conditions. In this case, boundary conditions (FF, PP, and FP) influence on natural frequencies was studied (as shown in Figs. 19,20). N m = 1 , 0 ≤ R ≤ 0.2m , x m,1 = x m = 0.35l , f j j = 1, 2, 3, . . . is the j-th natural frequency, and f j O is the j-th natural frequency when R = 0(J m = 0). Figure 19 shows that, regardless of the boundary conditions, each natural frequency becomes smaller as the radius of the concentrated inertia becomes larger, which suggests that the conclusions obtained in "Position  www.nature.com/scientificreports/ influence of the concentrated inertia", "Value influence of the concentrated inertia", "The concentrated inertia located at a certain point" and "Number influence of the concentrated inertias" sections are not affected by boundary conditions. In Fig. 20, it can be found that when J m is not located at x m = 0.35l , every mode shape, regardless of the boundary conditions, is changed with the increase of J m , the range of every modal is beyond the interval [− 1,1]. The findings means that the mode shape will be affected by the magnitude of the concentrated inertia.

Discussion of some phenomena
In "Application of AM" section, five interesting phenomena were observed. Phenomenon 1: When the concentrated inertia is located at x m = l n 1 (n 1 = 2, 4, 6, . . .) , some natural frequencies do not change with the magnitude of concentrated inertias; Phenomenon 2: As R(J m ) → ∞ , each natural frequency tends to a certain constant; Phenomenon 3: When the concentrated inertia is located at x m = l/n 1 (n 1 = 2, 4, 6, . . .) and R(J m ) → ∞ , the n 1 (2q + 1)/2 + 1-th natural frequency equals to the n 1 (2q + 1)/2-th natural frequency; Phenomenon 4: f max, k,n , f min, k,n in "Position influence of the concentrated inertia" and "Value influence of the concentrated inertia" sec-   Explanation of Phenomenon 2. Phenomenon 2 is caused by boundary conditions. As R(J m ) → ∞ , the position where the concentrated inertia is located becomes fixed boundary conditions. To explain the reason, Fig. 21 was given.
When R increases, the inertia of each concentrate inertias increases, and each of them is harder to be rotated. When R → ∞, then J m → ∞ , the position of each concentrate inertias is equivalent to the fixed end, and the shaft is divided into beam segments. Thus, the natural frequencies of the shaft mentioned in the title consist of the solution sets of all beam segments. Take the model in Fig. 21a for example. A shaft carrying single concentrated inertia J m located at x m = l/4 . When R → ∞ , J m → ∞ . As a result, the equivalent structure in Fig. 21b can be obtained. When the value of each parameter is set the same as in case three, Table 4 appears. From Table 4, the conclusion is proved and this conclusion is also applicable to multi-inertia systems. U 1 = ω 1 = n 1 π l n 1 = 0, 1, 2, . . . . U 2 = ω 2 = (2n 2 + 1)π 2x m |n 2 = 0, 1, 2, . . . . Explanation of Phenomenon 4. Maximum (minimum) values ( f max, k,n , f min, k,n ) in "Position influence of the concentrated inertia" and "Value influence of the concentrated inertia" sections appear alternately when the position of J m is changed. The position of each maximum value f max, k,n has no connection with J m , which is caused by the reason in "Explanation of Phenomenon 1" section. If ω 1 = ω 2 or ω 1 = ω 3 , the solution is x m = (2n 2 + 1)l/2n 1 or x m = [2(n 1 − n 3 ) − 1]l/2n 1 . Considering 0 ≤ x m ≤ l , n 1 > n 2 ≥ 0,n 1 > n 3 ≥ 0 and n 1 = 0 , the position of f max, k,n obtained was shown in Table 5, which is consistent with the phenomenon in "Position influence of the concentrated inertia" and "Value influence of the concentrated inertia" sections. However, the magnitude and position of f min, k,n changes with J m . Magnitude change of f min, k,n is caused by Magnitude change of J m , which had been given in "Explanation of Phenomenon 3" section. The position change of f min, k,n is caused by the combination of changes in the magnitude and position of J m . When J m → 0 and x m changes within 0 ≤ x m ≤ l , f min, k,n = f max, k−1,n . To clearly illustrate this conclusion, Figs. 22 and 23 was given. In this case, N m = 1 , N = k , and values of other parameters are set the same as in case three. Figure 22 was obtained by AM, which visually reflects the relationship between f min, k,n and f max, k−1,n . In order to prove the conclusion of Fig. 22, LMM was used for calculation ( R = 1 m ), and the results supported the conclusion in Fig. 22.  Table 5. Maximum position. (I) The torsional vibration model established in this study is universality. The free torsional vibration of shafts, such as ship shafting, rail car shafting, submarine propulsion system shafting, and other mechanical devices, can be simulated accurately by this model. (II) The analytical solutions of torsional vibration are given, which can be used in analyzing the natural frequencies and mode shapes of such structures and are instructive to the subsection design of shafting. The natural frequencies and modes vary with the magnitudes, positions, and the number of concentrated inertias. This provides a reference for the vibration reduction design of shafting. (III) It has been found that the results of AM are well in agreement with those of LMM, FEM, German, and Chen. The calculation method is very simple, which may be one of the simplest tools for studying the title problem.

Data availability
The data used to support the findings of this study are included in the article.